Scaling Relations for Watersheds 
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We study the morphology of watersheds in two and three dimensional systems subjected to differ- 
ent degrees of spatial correlations. The response of these objects to small, local perturbations is also 
investigated with extensive numerical simulations. We find the fractal dimension of the watersheds 
to generally decrease with the Hurst exponent, which quantifies the degree of spatial correlations. 
Moreover, in two dimensions, our results match the range of fractal dimensions 1.10 < df < 1.15 
observed for natural landscapes. We report that the watershed is strongly affected by local pertur- 
bations. For perturbed two and three dimensional systems, we observe a power-law scaling behavior 
for the distribution of areas (volumes) enclosed by the original and the displaced watershed, and 
for the distribution of distances between outlets. Finite-size effects are analyzed and the resulting 
scaling exponents are shown to depend significantly on the Hurst exponent. The intrinsic relation 
between watershed and invasion percolation, as well as relations between exponents conjectured in 
previous studies with two dimensional systems, are now confirmed by our results in three dimensions. 

PACS numbers: 64.60.ah, 91.10.Jf, 89.75.Da, 92.40.Cy 



I. INTRODUCTION 

Watersheds are the lines separating adjacent drainage 
basins (catchments). They play a fundamental role in 
water management [H-Q, landslides 043 1 ^^'^ flood pre- 
vention 0-11] • Natural watersheds are fractal [l^]- Ge- 
ographers and geomorphologists have found the evolu- 
tion of watersheds to be driven by local events classi- 
fied as stream captures or drainage rearrangements. These 
events can affect the biogeography [ll|, and may occur 
due to various mechanisms like erosion p^ - [T^ . natural 
damming Ml , tectonic motion [l5l - [l7| , as well as volcanic 
activity [18l |. 

Fractality of watersheds was first claimed in Ref. jTo| . 
Their observations were limited to small scales and only 
few samples. In Ref. extensive numerical simula- 

tions were devised to study watersheds on uncorrelated 
artificial landscapes, as well as on large-scale natural 
landscapes in the form of digital elevation maps (DEM) , 
as the ones obtained from satellite imagery [20] • Defin- 
ing the watershed to be the line dividing the entire land- 
scape into two parts, a novel and efficient identification 
algorithm was developed. Using this method, the self- 
similarity was confirmed and the fractal dimension was 
estimated to be 1.10 ± 0.01 for the Alps, 1.11 ± 0.01 for 
the Himalayas, and 1.211 ± 0.001 for uncorrelated arti- 
ficial landscapes. Due to the ubiquity of the obtained 
fractal dimension, also relations to other physical models 
have been proposed, such as optimal paths and optimal 
path cracks [21 23] , bridge percolation [24- 26,], and the 
surface of explosive percolation clusters [23, [2^. This 
opens a broad range of possible implications and appli- 
cations of the properties of watersheds. In this paper, we 
study the effect of correlations on the fractal dimension 
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of watersheds for artificial landscapes and compare in 
two dimensions to the values obtained for several natural 
landscapes ranging from fiat (Big Lakes) to mountainous 
regions (Rocky Mountains). 



The mechanisms triggering the stream capture and 
drainage rearrangement events 0, [T2| - [l8j . although di- 
verse, are all modifications of the topography. In that 
perspective, Fehr et al. [2§] investigated the effect of 
such events by applying small local perturbations to nat- 
ural and artificial landscapes and found that watersheds 
can be strongly affected. Power-law scaling behavior was 
observed, e.g. for the distribution of the area enclosed 
by the original and the displaced watershed, and shown 
to be independent on the strength of the perturbation. 
Tuning the correlation degree in the artificial landscapes, 
the values of the exponents were matched with the ones 
obtained for natural landscapes. Additionally, for un- 
correlated artificial landscapes, relations for the scaling 
laws to properties of invasion percolation j30| were con- 
jectured. We now extend these studies to three dimen- 
sional systems, where the watershed is a surface dividing 
the volume into two parts. We focus mainly on the study 
of conjectured scaling relations, analyze the finite-size ef- 
fects, and clarify the relation to invasion percolation. 



This paper is organized as follows. In Section |TT] we 
briefly revisit the deflnition of our model and show the 
dependence of the fractal dimension on the landscape cor- 
relations for both two and three dimensions. Section IlIII 
summarizes our results for the power-law distributions 
found for the effects of topographical modifications [2^ . 
To validate the relations conjectured for two dimensions, 
we study in Section IIVI the effect of perturbations for 
three dimensional systems. Conclusions are drawn in 
Section IVl 
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FIG. 1. (color online) Example of a watershed for uncorrelated systems with linear size L — 129 lattice sites in (a) two and 
(b) three dimensions. The upper and lower catchments drain to the top and bottom boundary, respectively. 



II. THE WATERSHED AND ITS FRACTAL 
STRUCTURE 

In nature, watersheds are the hnes separating adjacent 
drainage basins (catchments). We study the watershed 
Hnes in two dimensions as well as their extension to three 
dimensions. We use real and artificial systems in the form 
of Digital Elevation Maps (DEM) , giving the elevations 
(resistances) on a square (cubic) lattice. We define as the 
watershed, the line (surface) dividing the entire system 
into two sets of sites, catchments, as shown in Fig. [1] 
Each site in a catchment drains to the same boundary 
of a chosen pair of opposite boundaries (top-bottom in 
Fig. [TJ of the DEM. Hereby, the drainage of a site is 
determined by following the local slope and filling the 
valleys (local minima), until eventually reaching one of 
the two boundaries. For the determination of this line 
we use an iterative invasion percolation procedure (IIP), 
as introduced in Ref. [l9| . 

Recently, watersheds have numerically been shown to 
be self-similar objects [l^. These structures are typically 
characterized by their fractal dimension df, which is de- 
fined through the scaling of the mass M , corresponding 
to the number of sites or bonds in the watershed, with 
the linear system size L, 

M - L'^f . (1) 

In Ref. [Ti'l the fractal dimension of watersheds for un- 
correlated artificial landscapes has been estimated to be 
df — 1.211 ± 0.001, which differs from the ones obtained 
for landscapes from satellite imagery shown in Tab. HI 
This is expected as in natural landscapes typically long- 
range correlations are present [3l|. Therefore, we study 
here the dependence of the fractal dimension on correla- 
tions for several artificial landscapes. 

Spatial long-range correlated distributions can be ob- 
tained with fractional Brownian motion (fBm) j32l. [33j. 



TABLE I. Fractal dimension of watersheds for natural land- 
scapes obtained from satellite imagery [2^. We added here 
the values presented in Ref. [l9l ] for the Alps and the Hi- 
malayas for completeness. The fractal dimensions are ob- 
tained using the yardstick method. The error bars are of the 
order of 2%. 



Landscape 


df 


Alpes 


1.10 [19] 


Europe 


1.10 


Rocky Mountains 


1.11 


Himalayas 


1.11 [19] 


Kongo 


1.11 


Andes 


1.12 


Appalachians 


1.12 


Brazil 


1.12 


Germany 


1.14 


Big Lakes 


1.15 



Similarly to previous studies [231 2S. S^-lil]. we use the 



Fourier filtering method [23l. Iss'lssl- 37 1. which allows to 
control the nature and the strength of correlations. A 
detailed description of this method can be found, e.g., 
in Ref. [2^. In brief, the desired correlated distri- 
bution can be introduced by generating Fourier coeffi- 
cients in the reciprocal space according to a power-law 
spectral density. For each frequency in the reciprocal 
space, we calculate these Fourier coefficients through a 
random phase in the interval [0 ; 27r) and an amplitude 

s -2H-d 

\/k'f + ■ ■ ■ + k'jj , where fcj are the frequency in- 

dices of the discrete Fourier transform, d the spatial di- 
mension, and H the Hurst exponent. The inverse Fourier 
transform is applied to obtain the distribution in real 
space. Finally, we normalize the spatial domain distri- 
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bution in the range [0 : 1] to represent the correlated 
topology, characterized by the Hurst exponent H. Three 
different categories of fBm surfaces can be distinguished: 
0<H < 1/2, H = 1/2, and 1/2 < F < 1 [H]. The cor- 
relations between the increments are persistent (positive) 
for 1/2 < H < 1, meaning that sites with similar height 
tend to cluster together, leading to rather smooth sur- 
faces. The opposite is true for < H < 1/2, where the 
correlations of the increments are anti-persistent (nega- 
tive), resulting in surfaces that seem to oscillate more 
erratically. For H — 1/2, the classical Brownian mo- 
tion is recovered, where the increments are uncorrelated 
but the obtained heights are still correlated. The un- 
correlated distribution of heights is solely obtained for a 
constant spectral density, i.e. H = —d/2 {H = —1 and 
H = —3/2 in two and three dimensions, respectively). 

In Fig. [2] we plot the dependency of the fractal dimen- 
sion of the watershed on the Hurst exponent, measuring 
the degree of correlations. A monotonic decrease of d/ 
with H is observed, in line with what was observed for the 
optimal path crack in d = 2 Considering the typi- 

cal range of 0.3 < H < 0.5 for natural landscapes [3l| . 
our simulation results in two dimensions are in agreement 
with the values for the fractal dimension of watersheds 
in natural landscapes measured from satellite imagery, 
as listed in Tab. HI The lines in Fig. [5] show the fractal 
dimension obtained for the watershed of the Alps (solid 
line) and the one close to the Big Lakes (dotted line), 
characterizing the range of values obtained for natural 
landscapes. 

We now extend the concept of a watershed to a three 
dimensional system, in which the values at the sites no 
longer represent heights but for instance resistances. The 
watershed is now a surface, as shown in Fig. [T}), that di- 
vides the system into two parts. Similarly, each site in 
one part drains to the same boundary of a chosen pair of 
opposite boundaries (top-bottom in Fig. [T]) of the DEM. 
The drainage of a site is determined by following the low- 
est gradient in the resistances and filling the regions of 
local minima, until eventually reachin g on e of the two 
boundaries. Again, the HP procedure [1^ can be used 
to determine the watershed surface numerically. Alike 
the two dimensional case, we use Eq. ([Ij to estimate the 
fractal dimension. For uncorrelated three dimensional 
systems, we find df = 2.48 ± 0.02. In Fig. [21 we show 
that the fractal dimension of the watershed for three di- 
mensional systems also decreases continuously with the 
Hurst exponent H. 



III. IMPACT OF PERTURBATIONS IN TWO 
DIMENSIONS 

So far we have considered that the properties of the 
landscape, and consequently of the watershed, are static, 
i.e., do not change in time. However, landscapes might 
change due to several phenomena such as erosion, tec- 
tonic motion, and volcanic activity. Such changes in the 
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FIG. 2. (color online) The fractal dimension df of the wa- 
tershed as a function of the Hurst exponent H of the system 
in two (red stars) and three dimensions (green open circles) 
according to left and right hand axes, respectively. In two di- 
mensions, each point corresponds to an estimate of the fractal 
dimension by fitting the power law in Eq. ([T| to the water- 
shed masses for system sizes L — {5, 17, 65, 257, 1025, 4097}. 
For each system size, the mass was averaged over 10^ land- 
scape realizations. In three dimensions, the system sizes 
L = {5, 9, 17, 33, 65, 129, 257} were used with the same 
number of realizations. The lines show the fractal dimension 
of the watershed obtained for the Alps (solid) and close to the 
Big Lakes (dotted), characterizing, according to the left hand 
axes, the range of estimates for natural landscapes as sum- 
marized in Tab. in In the typical range of natural landscapes, 
0.3 < H < 0.5, the simulation results are in agreement with 
the natural ones. 



landscape are known to trigger local events called stream 
capture [l^ - [l8| , which can affect the watershed. In what 
follows, we extend the recent work presented in Ref. . 
This approach is based on a perturbation scheme, where 
a local event is induced by changing the height hk at a 
single site k of the system to hk + A, with A being the 
perturbation strength. It has been shown in Ref. [2^ 
that, without loss of generality, we can fix the perturba- 
tion strength to the height difference between the highest 
and lowest height of the landscape. As shown in Fig. [3]d, 
we quantify this response by the area, i.e. the number 
of sites, Ns, enclosed by the resulting watershed (dark 
(blue) line) after the perturbation at site k (cross) and 
the watershed of the original landscape (light (red) line) . 
The water can only escape from this area through one 
single site, which we call outlet. In this scheme, two out- 
lets can be considered, one before (cross) and another 
after the perturbation (dot). The former always coin- 
cides with the perturbed site k and can be connected to 
the latter, inside the enclosed area, by an invasion perco- 
lation cluster, whose mass we denote by M. The number 
of enclosed sites Ng, the mass M of the connecting clus- 
ter and the distance R between the two outlets (dotted 
line in Fig. |3^) are measured. After that, the original 
landscape is restored by resetting the height at k to its 




FIG. 3. (color online) (a) Example of a displacement (top, green) of the original watershed (bottom, red, the same as in 
Fig. [T^) after a perturbation at the spot marked by a cross for an uncorrelated artificial landscape. The dotted line connects 
the two outlets of the area before (cross) and after (dot) perturbation and defines the distance R. (b) The watershed close to 
the big lakes in the United States (light (red) line). A perturbation of 2 m at a spot (cross) near Thunder Bay caused a change 
in the watershed (dark (blue) line). The watershed displacement encloses an area of about 3730 km^ ;29]. The dot marks the 
new outlet of the area after the perturbation. 



initial value. Except of those sites located on the orig- 
inal watershed, this procedure is repeated for every site 
k in the landscape. In the following, we consider only 
perturbations actually leading to a displacement of the 
watershed, i.e. Ng > 0. To reduce finite-size effects, we 
explicitly exclude those perturbations, where the changed 
areas touch the borders of the system and, therefore, the 
original and the perturbed watershed are always overlap- 
ping at the boundaries. From the obtained set of mea- 
sures Ng, M, and R, we calculate the distribution P{Ns) 
of the numbers of enclosed sites (areas) Ng, the distri- 
bution P{M) of the clusters mass, and the probability 
distribution P{R) of the Euclidean distance R between 
the two outlets. To investigate the dependence of Ng and 
M on the distance i?, we define the average {Ng) and 
distribution P{Ng\R) of areas associated with a distance 
i?, as well as the average (M) and distribution P{M\R). 
For the dependence of M and Ng on each other we study 
the average mass associated with an area Ng. All these 
distributions and measures were sampled for each config- 
uration, then averaged over 2000 realizations of systems 
with size L = 513 and over 4000 for L = 129 and 257. 
In Ref. di], the distributions P{Ng), P{R), P{Ng\R), 
and P{M\R), as well as (Ng), have been shown to follow 
power laws of the form 



P{Ng)^N-^, 
P{R) - R-P, 

{Ng) ^ i^^ 

PiNg\R)^N-^, 
P{M\R) - Af~(i+"'\ 



(2a) 
(2b) 
(2c) 
(2d) 
(2e) 



and /? were found to increase with increasing H, the ex- 
ponent a did decrease. It has been observed that the av- 
erage angle between the lines connecting the outlets with 
the center of mass of the area, decreases with increasing 
correlation degree. This implies that, on average, the two 
outlets approach each other with increasing H, so that R 
is no longer representative of the area extension. There- 
fore, with 7? — >■ 1, fixing the distance R no longer re- 
stricts the areas entering the distribution P{Ng\R), such 
that the exponent a decreases with H and approaches 
p. For all considered landscapes the value a = 2 has 
been obtained. The distribution P{M\R) for uncorre- 
lated landscapes can be directly related to point-to-point 
invasion percolation, with a* = 1.39 being the subcritical 
exponent as introduced in Ref. [soj . 

In the following we present a functional description 
for the average number of enclosed sites {Ng) {R, L) that 
captures the dependence on i?, including also the finite- 
size effects. We can write the average for a given distance 
R as the first moment of the distribution P{Ng\R), i.e.. 



{Ng){R)^ / dNgNgP{Ng\R) 



where the exponents a, /3, and p have been estimated for 
several natural as well as artificial landscapes. While p 



(3) 



The distribution obtained from our simulation has a 
lower and an upper cutoff, as can be seen from Fig. |4^. 
Considering that those perturbations with areas touch- 
ing the boundary are excluded, we need to determine the 
scaling of the cutoffs numerically. From the data col- 
lapse achieved using the scaling P{Ng\R) = L^"/[A^sL^], 
we can see that the upper cutoff indeed scales as L'^. 
Similarly, the lower cutoff follows R^^ as obtained from 
the data collapse shown in the inset of Fig. |4^. Applying 
both cutoffs as bounds for the integral on the right hand 
side of Eq. ([3]), as well as for the normalization integral 
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1000 



FIG. 4. (color online) (a) Data collapse for the distribution 
P{Ns\R) with R = 1 (filled) and 7? = fO (open) and system 
sizes L — {f29, 257, 5f3} (triangles, circles, and squares, re- 
spectively), using the scaling P{Ns\R) = L'^" f[NsL'^]. The 
solid lines represent fits to the data of a power law with expo- 
nent a = 2.23±0.03. The inset shows the collapse of the same 
data for L — 513, when the x-axis is rescaled by R'', visual- 
izing the lower cutoff from the power-law behavior, (b) The 
average number of sites (area) (Ns) enclosed by the original 
and the perturbed watershed at an outlet distance R for the 
same landscapes and system sizes as in (a). The lines show 
the expression given by Eq. (Q. Each data point, in both (a) 
and (b), corresponds to an average over 2000 realizations of 
linear system size L = 513 and over 4000 for L — 129 and 
257. The error bars are smaller than the size of the symbols. 



for P{Ns\R), we obtain 



l-a 



(4) 



dA^., n: 



c 



l-a\ /L'i(2-") _i?d(2-«) 



This result matches, within the error bars, our sim- 
ulation data as shown in Fig. HJj. In the limit 
L ^ GO, only possible because a > 2 [2^, it reduces to 



{Ns) (R) = C(l - a)/ {2 - a)R'' and, therefore, ct = 2, as 
observed numerically. 

In a similar way, as calculated for the number of en- 
closed sites Ns , we investigate the mass AI of the invasion 
percolation cluster which connects the two outlets before 
and after perturbation. The cluster is always a subset of 
the enclosed sites Ns- Hence, we assume also the distri- 
bution P{M) and the average mass (M) to obey power 
laws of the form 



and 



P{M) - Af-(i+'5*) 



(5a) 
(5b) 



where we choose Eq. (j5ap to have a similar form as 
Eq. ((2e)) . Analogously to Eq. dS]), we can define the av- 
erage mass (M) as the first moment of the distribution 
P{M\R). With a (finite-size) scaling analysis similar to 
the one used to get the cutoffs in P{Ns\R), we find the 
lower and upper cutoff of P{M\R) to scale as R'^'f and 
L'^'f , respectively, where d*^ is the fractal dimension of 
the invasion percolation clusters. Using Eq. (l2e|) together 
with Eq. with M instead of Ng we obtain 



(Af) {R,L) = 



dMM- 



(6) 



dM M- 



-a*-l 



R f 



= c 



1 



j^d}{l-a') _ j^d}{l-a') 



In Fig. [5] the matching, within the error bars, of this 
result with the simulation data is shown. For L — oo the 
resuh reduces to (M) (i?) = a* /{a* - l)R'^f , as a* > 1, 
yielding cr* = = 91/48. The relation of the area 
and the mass is analyzed with the average mass (M) 
associated to a given number of enclosed sites Ng. From 
the dependence of (Ng) and (M) on R we expect 



(7) 



which matches the value 0.95 ± 0.01 obtained from the 
power-law fit to the data plotted in the inset of Fig. [5] 
Based on this result, using Eq. ([7]) together with Eq. (I^al) 



and (Hi]) we expect /3*=(T/3/dJ-lsa0.22. This predic- 
tion matches the value (3* ~ 0.21 ± 0.02 estimated from 
the data collapse of P{M) shown in Fig.[6l Summarizing 
our results in two dimensions, together with the relations 
conjectured in Ref. |29j we obtain the relations 



a 



zL 



/3=^(l + /5*), 
P = crl3, 



(8a) 

(8b) 
(8c) 



between the exponents a and f3 for the numbers of sites 
Ns enclosed by the original and the perturbed watershed. 
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FIG. 5. (color online) The average mass (M) of the inva- 
sion percolation cluster connecting the two outlets with a dis- 
tance R for uncorrelated two dimensional landscapes of sizes 
L = {129, 257, 513} (triangles, circles, and squares, respec- 
tively). The lines show the predictions according to Eq. © for 
different system sizes. The inset shows the average mass vary- 
ing with the number of enclosed sites. The line is a power-law 
fit to the data yielding an exponent 0.95 ±0.01 ~ d*f/2. Each 
data point corresponds to an average over 2000 realizations 
of systems with size L — 513 and over 4000 for L = 129 and 
257. The error bars are smaller than the size of the symbols. 




* 



M/L"/ 

FIG. 6. (color online) Data collapse of the size distribution 
P{M) of the mass A'l of the invasion percolation cluster con- 
necting the two outlets for uncorrelated landscapes of three 
different system sizes L = {129, 257, 513} (triangles, circles, 
and squares, respectively). The line represents a power-law 
fit to the data for the largest landscape (squares) yielding an 
exponent 1 + P* — 1.21 ± 0.02. Each data point is an average 
over 2000 realizations of systems with size L — 513 and over 
4000 for L — 129 and 257. The error bars are smaller than 
the size of the symbols. 



a* and /3* for the mass M of the invasion percolation 
clusters, and p for the probability density to induce, af- 
ter perturbation, a change at a given distance. Despite 
the similarity of a and p for uncorrelated two dimen- 
sional landscapes they differ significantly in three 



dimensions, as shown in the next section. 



IV. IMPACT OF PERTURBATIONS IN THREE 
DIMENSIONS 

In this section we extend the previous concepts to three 
dimensional systems. As introduced in Sec. the wa- 
tershed in three dimensions is a surface that divides the 
system into two parts. Similar to two dimensions, a per- 
turbation is induced by changing the local resistance 
at site fc to r/c -f A, where A is the perturbation strength. 
We quantify the impact on the watershed surface, as be- 
fore, by the number of sites Ng enclosed by the original 
and the perturbed watershed, which now corresponds to 
a change in volume. As an example, the largest of all 
changed volumes observed in our simulations is shown 
in Fig. [71 The outlets of this volume, before and after 
perturbation, are determined as well as their distance R 
and mass M of the invasion percolation cluster connect- 
ing them. The original system is then restored by reset- 
ting rk to its initial value. This procedure is repeated 
for all sites k of the system, except those located on the 
watershed. Again, the averages (Ng), (M), and distri- 
butions P{R), P{N,), P{N,\R), P{M), and P{M\R) are 
sampled and averaged over 2000 configurations. Similar 
to the two dimensional case, as shown in Fig. [5] for the 
distribution P{Ns) of volumes, we find these quantities 
to follow power laws of the form introduced in Eqs. ([2a|- 
(Pe)) and Eqs. (|5ap - (j5b[) . The values of all the exponents 
estimated for uncorrelated systems are summarized in 
Tab.ni 

TABLE II. Summary of the exponents numerically obtained 
in this study. Estimates are given for watersheds on natural 
and uncorrelated artificial landscapes in two dimensions, as 
well as for uncorrelated artificial systems in three dimensions. 
For exponents obtained in previous works the corresponding 
citations are given. The range 1.10— 1.15 of fractal dimension 
for natural landscapes is obtained from Tab. IH 



Type 


Eq. 


natural 


artificial (uncorr.) 


d 




2 




2 




3 


df 


O 


1.10 - 1.15 




1.211±0.001 [19] 


2.48±0.02 


a 


(|2d|| 


2.3±0.2 


[29] 


2.23±0.03 




2.4±0.1 


P 


(|2a|| 


1.65±0.15 


m 


1.16±0.03 


[29] 


1.31±0.05 


P 


(|2b| 


3.1±0.3 


[29] 


2.21±0.01 


[29] 


3.2±0.2 


a 


((2^ 


2 


[29J 


2 


[29] 


2.45±0.05 


a* 


(I2il 






1.39±0.03 


[29, 30] 


1.4±0.1 




(|5a|| 






0.21±0.02 




0.29±0.06 


a* 


(|5b| 






91/48 


[46] 


2.53 [46] 



Although the volumes Ng are still compact, we find 
cr = 2.48 ± 0.02 by an analysis of (Ng) and P(7Vs|i?) 
similar to the one performed in two dimensions using 
Eq. ^ in the limit L oo. The obtained value of a 
is close to the fractal dimension of the watershed, what 
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FIG. 7. (color online) (a) The largest of all observed changed volumes (dark, blue) attached to the perturbed watershed (light, 
red) for the same uncorrelated three dimensional system of linear size L = 129, as used in Fig. [TJd. The upper and lower 
catchments drain to the top and bottom borders, respectively, (b) The same changed volume as in (a) without the watershed 
for better visibility. The lines mark the sites where the original watershed intersects the system boundaries. 




FIG. 8. (color online) Data collapse of the distribution P{Ns) 
of the number of sites (volume) enclosed by the original and 
the perturbed watershed for uncorrelated three dimensional 
systems of size L — {33, 65, 129}, (squares, circles, triangles, 
respectively). The line represents a power-law fit to the data 
yielding an exponent /3 = 1.31 ± 0.05. Each data point is 
an average over 2000 realizations. The error bars are smaller 
than the size of the symbols. 



suggests that the mass of the changed volume is domi- 
nated by its surface, in contrast to what we observe in 
two dimensions. This is also confirmed by a box counting 
analysis of the largest of all changed volumes, as depicted 
in Fig. [5] In agreement with the findings by Lee in Ref. 
p7| . that the size distribution of sub-critical point-to- 
point invasion percolation, as introduced by Araiijo et 
al. [30| , is independent on the dimensionality of the sys- 
tem, we obtain a* = 1.4 ± 0.1 for uncorrelated systems, 
similar to the value found in two dimensions. This con- 



FIG. 9. (color online) Number of cubic boxes C of size e 
covering the changed volume, shown in Fig. [71 as obtained 
from a box counting method. Each data point consists of a 
single measurement. The line corresponds to e^ **. 



firms the relation to invasion percolation drawn in Ref. 
[2£|. Inserting the estimates of dj, a, and a* , given in 
Tab. [m into Eq. dHi]), we find a w 2.48 matching with 
the value 2.4±0.1 obtained numerically, hence validating 
the conjectured relation. Similarly, we also observe that 
our results are consistent with Eq. (I8bp . 

The dependence of the exponents a, /?, and p on the 
Hurst exponent H is shown in Fig. [TUl confirming the 
relation p = a/3 (Eq. (jSc]) ). independent on the degree of 
correlations. We observe a to decrease for H > and to 
approach /?, similar as found in two dimensions (compare 
Section mil. 
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FIG. 10. (color online) The exponents a, j3, p, a (squares, 
circles, triangles up, and triangles down, respectively) as a 
function of the Hurst exponent H for perturbations in 3D. 
Each data point consists of a similar measurement as per- 
formed in Fig. [5] to obtain /? for uncorrelated systems. 



CONCLUSION 



In summary, we were able to show that watersheds 
are fractals with a fractal dimension in the range of 
1.10 < df < 1.15 for all analyzed natural landscapes 
from mountainous (e.g. Alps) to rather flat ones (e.g. Big 
Lakes). By studying model landscapes with long-range 
correlations characterized by the Hurst exponent H, we 



determined the dependence of the fractal dimension on 
H in both two and three dimensions, where for the for- 
mer we found good quantitative agreement with natural 
landscapes, for which 0.3 < H < 0.5. Extending the work 
done in Ref. ^29,] on the impact of perturbations on wa- 
tersheds in two dimensional systems, we found the lower 
and upper bounds of the distribution of areas to scale 
with distance and system size with an exponent cr = 2 in 
two and a = 2.48±0.02 in three dimensions. Considering 
these cutoffs, we were able to derive a function describ- 
ing the average area, which includes the finite-size effects. 
In a similar way, by integrating P(M\R), we obtained a 
function for the average mass (M) of the invasion perco- 
lation cluster. In the extension of the perturbation study 
to three dimensional systems, we observed the changed 
volumes to be dominated by their surface, i.e., the wa- 
tershed, resulting in a value of cr close to the fractal di- 
mension of the watershed itself. The distributions and 
averages for the three dimensional case are shown to fol- 
low power laws as in two dimensions. Indeed, we found 
in the uncorrelated case an intrinsic relation to invasion 
percolation. Finally, our results in three dimensions are 
consistent with the conjectured relations between the ex- 
ponents. As a foUowup of the work presented here, it 
would be interesting to see how the impact of perturba- 
tions could be related to other physical models, such as 
optimal path cracks [2l|,[2l], bridge percolation ^3, and 
the surface of explosive percolation clusters [27|,128|. 
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and thank CNPq, CAPES, FUNCAP, and the 
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